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Abstract: This study proposes a fast 3D dynamic programming expansion to find a 
shortest surface in a 3D matrix. This algorithm can detect boundaries in an image 
sequence. Using phantom image studies with added uniform distributed noise from 
different SNRs, the unsigned error of this proposed method is investigated. Comparing the 
automated results to the gold standard, the best averaged relative unsigned error of the 
proposed method is 0.77% (SNR = 20 dB), and its corresponding parameter values are 
reported. We further apply this method to detect the boundary of the real superficial femoral 
artery (SFA) in MRI sequences without a contrast injection. The manual tracings on the 
SFA boundaries are performed by well-trained experts to be the gold standard. The 
comparisons between the manual tracings and automated results are made on 16 MRI 
sequences (800 total images). The average unsigned error rate is 2.4% (SD = 2.0%). The 
results demonstrate that the proposed method can perform qualitatively better than the 2D 
dynamic programming for vessel boundary detection on MRI sequences. 
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1. Introduction 

Many image edge detection methods have been proposed in the last two decades [1,2]. They mainly 
differ in the types of smoothing filters that are applied and how the edge strength measures are 
computed. Edge detection methods allow broken or discontinuous edge lines. Boundary detections, 
however, are somehow different, as some applications do not allow discontinuous boundaries. Many 
types of methods can handle the vessel extraction problem. Several hundred papers are surveyed in [3], 
which consider vessel extraction in two-dimensional images. For instance, two-dimensional boundary 
detection utilizing graph-searching principles [4,5] has been studied frequently in medical image 
analysis. Boundary detections are sometimes transferred to an optimization problem that minimizes a 
given cost function [6-9]. Dynamic programming (DP) is an optimization tool often used in boundary 
detection [6,8,9]. Dynamic programming method performs well if the boundaries are visible. The 
motivation of this study is raised from the observation that the femoral artery in the cross-sectional 
view in MRI sequences without contrast injection is sometimes invisible. The artery is visible if the 
blood flow velocity is large enough that it appears in contrast compared to the surrounding tissue. This 
is most often an MRA image sequence. If the blood flow velocity is small during a short time period, 
the artery has limited contrast and is not visible. Under this extreme situation, the 2D DP fails to detect 
its boundary because it is hard to obtain a feature, normally the gray-level gradient, to represent the 
boundary. Our previous study has applied the local contrast to add additional information to handle 
extreme cases [10]. It still used local information but not information between two succeeding images. 
A relative study proposed a graph-based method to find an optimal surface in volumetric images [11]. 
This article is most related to our method, but it differs in principal. Good reviews are reported in other 
studies [12,13], in which methods are considered or models are represented in three dimensions. The 
methods treated the three-dimensional vessel geometry problems, in their issues vessel branching is a 
major problem to be solved. This is, however, not the major problem in our study. We do not consider 
three-dimensional vessel geometry. Instead, we focus on the vessel cross-sectional lumen changing with 
respect to a heart cycle at the same place. The lumen area must be quantified accurately. Moreover, 
the MRA images without a contrast medium sometimes present difficulty when visualizing vessel 
boundaries. Previous methods [3,12,13] are not suitable to solve our problem. This is for two reasons: 
(1) Our MRI data have vessel images, in which some of them the boundaries are vague and hard to 
recognize, even by human beings; (2) Our goal is to achieve a high accuracy on vessel's crossing area 
quantization. The goals of previous 3D methods are more suitable for visualization on vessels 3D 
structure, although the quantification is also possible if vessel boundaries are obtained. 

To solve this problem, we propose a 3D-expansion of DP that can seek an optimal surface in a 3D 
matrix. This method can search succeeding boundaries in an image sequence, as succeeding boundaries 
form a surface in a 3D matrix. The relationship between two succeeding boundaries is thus considered in 
the smoothness constraint. To investigate the performance of this method, we design experiments to 
calculate its accuracy. The significance of this paper is two-fold: (1) we propose a 3D-expansion of DP 
which can locate an optimal surface among a 3D matrix; (2) we explore the best parameter set of the 
proposed 3D expansion of DP for round shape objects' boundary detection, and its unsigned relative 
error is reported. 
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The rest of this paper is organized as follows. In Sections 2 and 3, we introduce the 2D and 
3D-expansion of DP, respectively. In Section 4, we address our experimental designs. In Section 5, we 
illustrate the experimental results. We then discuss the proposed method and give a conclusion in 
Sections 6 and 7. 

2. Methods 

2.1. 2D Dynamic Programming 

In computer science, DP is a method for solving complex problems, which has similar properties 
that can be broken down into simpler sub-problems. Therefore, to solve a problem is then transformed 
to solve different parts of its sub-problems. However, many sub-problems are actually the same. The 
DP solves each sub-problem once and the result is saved for solving next similar sub-problem. Thus 
DP reduces the number of computations. A good tutorial on DP can be found in [14]. This approach 
was often used to find the contour of an object such as a road [15], and it has many applications on 
boundary detection in medical images [6,8-10,16-18]. A study has shown that DP outperforms some 
edge detectors against noise, including the Canny and Prewitt detector [19]. If the object to be detected 
contains noise, the DP might therefore be a good detection method, as for an MRI breast mass 
segmentation [20]. Similarly, DP is suitable for boundary detection in noisy images, including 
sonography [8,21] and MRI vessel segmentation without a contrast injection [10]. The 2D boundary 
detection using DP can be transformed to an optimization problem seeking an optimal path in a feature 
image. The searching direction is performed from one side to its opposite side. Converting a closed 
contour from its center in the Cartesian coordinate to the polar coordinate is thus required [22] for 
closed contour detection. 

Let the feature image saved in a 2D matrix be F, F e 9t Mx7V . The boundary detection problem is 
then transformed to an optimization problem that searches for an optimal path. Assume the searching 
direction is from left to right. The optimization function can be defined to find the optimal path having 
the following global minimum: 

7>argmin^ =i {F(j;J|x = l,2,3,---7V} (1) 

subject to some constraints, where y x is the y-coordinate on the x-th column in matrix F, and j/ x and j/ x +i 
are neighborhood. This optimization function can thus be reformulated to implement DP with respect 
to an iterative cost function: 

C(x, y) = min M _ dd) C(x -l,y + j) + F(x, y) + a\j\ (2) 

subject to 1 < x < N 9 1 < y < M, where a is a weighting parameter controlling the smoothness of the 
searched path and d is the maximum distance between two connected nodes. C{xy) is a two-dimensional 
cost map. The global optimization problem is the same as its subproblem C{x - ly\ C(x - 2y) and 
vice versa. We set C(ly) = F(ly\ as a boundary condition. The optimal index f can be determined by 
the following equation: 

/ = arg min ye ( _ d d) C(x -l,y + j) + a\j\ (3) 
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The index can be stored in the 2D coordinate matrix Y(xy) —y +j* 9 which is a pointer indicating a 
point on the previous column (x - 1). The cost map and path links are thus constructed column- wise 
from left to right on the feature matrix F. After construction, the optimal path can be found by tracing 
the path link backwards on the last column (x = N), which has the global minimum. There are some 
notable variations of the 2D DP, including the dual- [8] and multi-path DP [23]. 

2.2. 3D-Expansion of Dynamic Programming 

Let the 3D matrix R have size M x N x P, where M and N are the numbers of rows and columns, 
and P denotes its depth. The volumetric matrix contains feature image sequences of interest. Values in 
R are normalized, i.e., 0 < i?(y,x,z) < 1, where y, x, and z are indices of the corresponding dimensions. 
Assume that features are saved in the 3D matrix and the goal is to find an optimal surface having the 
shortest path and lowest value summation from one side on the z-axis to another side with some given 
constraints. Figure 1 shows the 3D matrix R and the surface to be detected. The typical constraint in 
DP controls the smoothness or continuity of the sought surface. Assume the searching direction is from 
x = 1 to x = N and z = 1 to z = P. Two parameters control the smoothness: d\ > \y x -y x -i\ controls the 
smoothness on the x - y plane TV x M, and di > \y z -y z -\\ controls the smoothness on the x-z plane Nx P. 
The parameters allow the maximum distance between two connected nodes. The ^-coordinate (y{x,z)) 
denotes the height of the optimal surface, which can be found via an optimization problem: 



subject to the constraints \y x -y x -\ \ < d\ and \y z -y z -\\ < d^. 

Figure 1. The 3D matrix R contains the feature of an image sequence. The structure of the 
accumulation matrices (a) Ci and (b) C2 is organized in different orientations. 



N P 



yl = arg min £ X )\^ x < N,l< z < P} 



(4) 



x=\ z=l 



Ci(y,x,z) 



c 2 (y,z, X ) 



y 




(a) 



(b) 



This optimization problem minimizes the iterative cost functions as such: 

C x (y, x, z) = min (Q (y + i 9 x — 1, z) + a x \ i |) + R( y, x, z) 



(5a) 



for 1 <x<7Vand 1 <z<P. 
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C 2 (y, z, x) = min (C 2 (y + i,z- 1, x) + a 2 \ i |) + C, (y, x, z) r 

-d 1 <i<d 1 V^ u / 

for 1 < z < P and 1 < x < TV, where CI and C2 are accumulation matrices of size M x N x P and 
Mx P x N, respectively. Figure 2 shows their organizations. The initializations of CI and C2 are thus: 

Ci(*,l,z) = R(*,l,z)for 1 <z<P (6a) 

C 2 (*, 1 ,x) = R(*, 1 ,x) for 1 < x < TV (6b) 

Figure 2. The accumulation matrices Ci and C2. The dimensions of (a) Ci and (b) C2 is 
same as those of R. 




(a) (b) 

During the accumulation process, the indices must be stored in Y\{x,y,z) for Ci and Y 2 ix,y,z) for 
C2 for the backwards tracing [8]. After the minimization process, the accumulation matrix C2 contains 
the minimal value in the last column. In the backwards tracing, the y-coordinates of the optimal 
surface can be obtained sequentially. The minimal value in the last column is found using 
y * (x, end) = min C 2 ( y, end, x) . 

y 

The remaining y- coordinates in the same index x can be traced backwards: y*(x,end - 1) = 
Y2(y*(x,end)jc 9 end). The node in Ci searches for a node in its previous column (x z - 1), which makes 
the accumulation minimal and saves the optimal path link in another 3D matrix Yi, shown in 
Figure 2(a). During this process, the connection is considered only at the same depth z. The possible 
searching range is denoted by d\ and weighted by a\\ both control the smoothness of the surface on 
xy-plane. Though the searching process is determined by the node value on its previous column plus 
the smoothness a\\i\ 9 as in Equation (2a), the final value of the node is summed by the node value in 
the feature matrix i?(y,x,z). Similarly, the node in C2 searches for the node in its previous column 
(zk - 1), which makes the accumulation minimal and saves the optimal path link in another 3D matrix 
Y 2 , shown in Figure 2(b). During this process, connection is considered only the same depth x. The 
possible searching range is denoted by ^2 and weighted by <i2\ both control the smoothness of the 
surface on j/z-plane. Though the searching process is determined by the node value on its previous 
column plus the smoothness a?\i\ 9 as in Equation (2b), the final value of the node is summed by the 
node value in Ci(y,x,z), in which the neighboring information on xj/-plane has been already considered. 
The final C2 matrices thus contain the feature and neighboring information on both xy- and j/z-planes. 
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Compared to the traditional DP, the 3D DP has the advantage that it connects curves in two planes 
(xy- and j/z-planes). This property allows extracting the boundary even if the boundary is vague. 
Section 5 presents the examples. Notably, if the backwards tracing process is performed on the Ci 
matrix, it produces the same results as the traditional DP. 

3. Experimental Design 

3.1. Phantom Study 

To investigate the accuracy of the proposed method, we design a phantom image sequence with 
added uniform distributed noise. The image gray level is normalized to be within the range [0,1]. The 
foreground (artery lumen) is set to 1, and the background is set to 0.4. The image has the size 41 x 41, 
with simulated artery lumens having radii ranging from 9 to 12 pixels (with step 1) and back to 8 pixels. 
Each image contains one artery lumen, and there are eight images in a sequence. The noise intensity is 
randomly given at every pixel, and the SNR (signal-to-noise ratio) is measured as such: 

SNR = -J— V V 20 log 10 (6S 

where M = 41 is the image size and g(x,y) and n(x,j/) represent the raw image gray-level and added 
noise (zero mean, uniform distributed) intensity at coordinate (x,j/), respectively, where |w(jcj;)| - n ievei- 
The SNR is controlled by ni eve h where 0 < ni eve i < 1 . The images are smoothed by a Gaussian function 
having a = 0.5 before noises are added, which makes images more realistic. 

To extract boundary information, we apply the directional gradient [6]. We investigate the optimal 
parameter, which is set to achieve the best accuracy of the proposed method for round shape boundary 
detection. There are two parameters to be explored: the image resize and weighting factors s = a\ = 
in Equations (2a) and (2b). We fix the constraints as d\ = 1 and d2 = 2. The factor d^ is larger because we 
allow the method to catch larger radius changes in the sequence. The image resize is performed after the 
image gradient computation, and it is transformed from Cartesian to polar coordinates. A 3D volume 
containing the gradient images is then used as input for the 3D-expansion of DP. 

3.2. MRA Real Images 

A mobile MRI (Siemens-type "AvantoTM", 1.5 T) was used for image acquisition [9,10]. Twelve 
of the forty-four participants in the TEFR09 study participated in this study after approval from the 
local ethics committee of Ulm University in accordance with the Declaration of Helsinki. Complete 
MRI sequences were acquired from all subjects. To validate the novel SFA (superficial femoral artery) 
lumen detection algorithm presented in this study, several MRI sequences from selected subjects were 
randomly selected. 

A mobile 1.5-T Magnetom (Siemens-Avanto™, Model Mob. MRI 02.05, Siemens Ltd., Erlangen, 
Germany), with a flexible 6-channel body matrix coil with 6 integrated low-noise preamplifiers 
(Siemens Ltd.) and bilateral table fixation, was used to acquire SFA MRI sequences from the subjects. 
All of the participants were athletes. The athletes were fixed in a stretched supine position, head 
forward on the MR table. To identify the axial perpendicular acquisition location at the right SFA 



Sensors 2012, 12 



5201 



10 mm beneath the bifurcation of the common femoral artery, a biplane coronal and sagittal localizer 
(TRUFI: "true fast imaging with steady state precision"; Siemens Ltd.) was used. As in the common 
carotid artery MR-measurement, changes in the SFA diameter during systole and diastole was assessed 
with a T2*-weighted gradient- spoiled gradient-echo cine-sequence (FLASH: "fast low angle shot", 
Siemens Ltd.) in a two-dimensional cross section view with prospective two-dimensional ECG gating 
(cardiac triggering). Specific sequence parameters were set thus: flip angle 15°, echo time variable 
between 4-6 ms (depending on heart rate), repetition time variable 20-40 ms (depending on heart 
rate), slice thickness 6 mm, field of view 768 cm 2 , matrix size 512 x 384, pixel width 0.625 mm (pixel 
area 0.3906 mm 2 ) ISO, pixel bandwidth 250, and 50 images per sequence for one RR-cycle 
(approximately 300 heart beats per sequence). The total imaging acquisition time was approximately 
4 min 30 s to 5 min for each sequence. Figure 3 shows two typical MRA images: 3(a) was taken in the 
systolic phase and 3(b) in the diastolic phase. 

Figure 3. Two raw MRA images in one sequence. Each sequence contains 50 images, 
(a) Image taken in the systolic phase; (b) Image taken in the diastolic phase. The arrow 
indicates where the SFA is. 





Each real MRA image sequence contains 50 images. Our previous study has proposed an automatic 
SFA position detection method [10]. An ROI (region of interest) can be extracted using that technique, 
in which the right SFA is in the center. The directional gradient [9] is then applied to the ROI to extract 
the boundary information. The extracted gradients are transformed from Cartesian to polar coordinates. 
Fifty total gradient images form a 3D volume to be input for the 3D DP. 

3.3. Accuracy Analysis 

The proposed system is applied, and the SFA cross-sectional lumen area of each image is calculated 
for the following comparison. The comparison is performed by calculating their relative unsigned 
errors: 

Si = | AManual(i) ~ A Automate JJ)\' 'A Manuali}) X 100% (8) 

where A Automate Jj) and A Ma nuai(i) are areas calculated by the automated and manual drawing on image 
number /, respectively. The mean errors (on 50 images) and standard deviations (SD) can be calculated. 
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4. Results 

Figure 4 illustrates phantom images with different noise levels (SNR). To investigate the effects that 
these two parameters (resize factor and weighting factor s) have on accuracy, four different phantom 
image sequences, with SNR = 14, 16, 18, and 20 dB, are tested. Figure 5 shows the corresponding 
averaged accuracies. The z-axis in Figure 5 denotes the averaged unsigned error percentage (for eight 
images and ten experimental repetitions). The value of s is discretized from 0.01 to 0.21 with a step of 
0.01. The resize factor value is discretized from 1 to 3 with a step of 0.1. From the results, we found that 
parameter s does not affect accuracy too much but the image-resize factor does impact accuracy. The 
optimal image resize factor value is 1.6, no matter the noise level. When SNR =14 dB, the unsigned 
error ranges from 1.50% to 2.47% with resize factor = 1.6. The best accuracy appears when s = 0.01. 
Without image resize (factor =1), the unsigned error ranges from 3.15%> to 4.2 1%>. The best accuracy 
appears when s = 0.02. When SNR = 20 dB, the unsigned error ranges from 0.11% to 1.62% with resize 
factor = 1.6. Without image resize (factor =1), the unsigned error ranges from 3.05%> to 3.46%>. The best 
accuracy appears when s = 0.06. Figure 6 demonstrates the contour detection results (SNR = 14 dB). 

Figure 4. The phantom images with added noise from different SNRs. 




(e) <0 (9) (h) 

Furthermore, we demonstrate the ability of the proposed method when one image is ruined. 
Figures 7 and 8 use the same eight phantom images (SNR =14 dB), and their boundary detection 
results differ on the ruined image. The traditional DP has no knowledge about the third direction, and it 
fails to detect the boundary if the image is ruined. However, our method shows its robustness against 
strong noise. 
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Figure 5. The unsigned error plots of different phantom images. SNR is (a) 14 dB, 
(b) 16 dB, (c) 18 dB, and (d) 20 dB. 




0 1 . 0 1 

resize factor b resize factor 



(c) (d) 



Figure 6. The contour detection results. SNR = 14 dB, s = 0.01, resize factor = 1.6. 
Averaged unsigned error = 1.5%. 
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Figure 7. One image is ruined. The proposed method can still hold the contour because of 
the continuity of the 3D DP method {di =1). 




In the real MRA image study, we show the robustness of the proposed method against vagueness. 
Figure 9 shows the results of nine sequential images and their corresponding detected boundaries. They 
are organized in a 3-by-3 matrix. In each element of the matrix, the upper image is the raw sub-image, 
and the lower image has its result superimposed on it. From the results, we can see that the artery's 
boundaries are vague in the first five images. It is even difficult for a medical expert to define their 
boundaries. Our algorithm can overcome this difficulty by using the 3D dynamic programming that 
has considered the connectivity information in three dimensions in these vague images. The proposed 
method thus has superior performance against vagueness on the boundary location. 
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Figure 9. Boundary detection results. Nine sequential resultant images are shown. The 
vessel boundaries are vague in the first five images. Their exact boundaries are difficult to 
determine. Based on the continuity consideration in three dimensions, the proposed method 
can detect the correct vessel boundaries. 
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The Bland- Altman plot [24] is often used in a clinical comparison of a new measurement technique 
with an established one to see whether they agree sufficiently for the new technique to replace the old. 
Figure 10 shows the comparison results between the proposed method and the experts' manual tracings 
via the Bland- Altman plot. Figure 10(a) is the result of one sequence having the mean error rate 

2 2 

2.1% ± 2.6%. The mean of the difference (y-axis) is -0.19 mm with SD 2.14 mm . The 95% confidence 
intervals are within 4.2 mm . Figure 10(b) is the result of 16 sequences (800 images) having the mean 

2 2 

error rate 2.4% ± 2.0%. The mean of the difference (y-axis) is 0.22 mm with SD 2.11 mm . The 95% 
confidence intervals are within 4.1 mm 2 . The bias is very limited (0.22 mm 2 ). The system is consistent 
in its data number (SD = 2.11 mm 2 ). The results demonstrate the system's ability to replace the 
expert's manual work. 
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Figure 10. Bland-Altman plots. The solid line indicates the mean, and dash lines are 
± 1.96 SD. (a) The plot of sequence no. 9 (50 images). (1.96 SD = 4.2 mm 2 ); (b) The plot 
of all 16 sequences (800 images). (1.96 SD = 4.1 mm 2 ). 
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Figure 11. The cross-sectional SFA area changing in time (sequence no. 9). The unit in the 
y-axis is mm 2 . The x-axis unit is image number. 

SFA cross-sectional area. 
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Figure 11 illustrates the result of sequence no. 9. Notably, the system can detect the artery's area 
decreasing (around image number 5 in this sequence), which is a physiological phenomenon. This is 
the most difficult issue, as the artery's boundary in this phase is unclear and it has the least image 
contrast. Sometimes even experts cannot exactly locate the boundary. Table 1 provides the accuracy 
analysis of the system using the unsigned errors rate of all 16 sequences. The maximum mean error 
rate is under 3.6%, and all SD are under 2.8%. 



Table 1. The average unsigned relative errors on 16 MRI sequences (800 total images). 



Sequence No. 


Mean (%) 


SD (%) 


Sequence No. 


Mean (%) 


SD (%) 


1 


2.4 


1.8 


9 


2.1 


2.6 


2 


3.2 


2.1 


10 


2.9 


2.5 


3 


1.4 


1.3 


11 


1.4 


1.0 


4 


3.2 


2.7 


12 


2.5 


1.7 


5 


3.6 


2.6 


13 


2.4 


1.2 


6 


3.1 


2.4 


14 


3.2 


2.0 


7 


1.9 


1.6 


15 


1.6 


1.7 


8 


2.7 


2.8 


16 


1.4 


1.1 


Mean 


2.4 




SD 


2.0 





Each sequence contains 50 images. Resize factor = 1.6, s = 0.1. 



The computer system has an Inter Core™ 2 CPU T5600, 1.83 GHz, with 2 GB RAM. All programs 
are designed based on the Matlab platform [25]. The computation time for each sequence (50 images) is 
approximately 2.1 s. The computation is much faster than our previous study [10], in which it takes 
approximately 25 s under the same hardware construction. 

5. Discussions 

In the proposed algorithm, there are two weighting factors, a\ and «2, and two constraint factors, d\ 
and d,2 (Equation (5)). The constraint factors allow the maximum neighboring boundary, whereas the 
weighting factors affect boundary continuity. Distance d\ defines the boundary connectivity on the 
same image. It is thus usually defined to be less than or equal to 2. Distance defines the boundary 
connectivity on a different image level and is thus case dependent. If boundaries on different levels 
have a large distance, then <i 2 must be defined larger. In this study, the SFA have 50 images during an 
R-R interval. The boundary thus does not change much, and di is set to 2. 

From the phantom study, we know that the effect of the resize factor is larger than the weighting 
factors a\ and «2- We therefore set the resize factor to 1.6 (the optimal parameter value in the phantom 
study) and choose the weighting factors as0.1(ai = a2 = s = 0.1) arbitrarily in the real MRI study. 

We have proposed a three dimensional DP, which can find an optimal surface in a 3D matrix where 
features are normally saved. This method has robustness against strong noise, especially when one or 
two images are absent in the image sequence. The proposed property is useful, especially for detecting 
an artery's boundary in an MRI when one or two images in the artery's boundary are vague. The 
connectivity can help rebuild the artery's boundary. This method can be extended to track an object's 
boundary in image sequences, especially when partial occlusion exists. Furthermore, it is possible to 
embed another model into the DP and let it extract a special shape, including a round [6] or oval shape. 
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The limitation of this proposed system is that it cannot handle branches. Many vessel segmentation 
methods must face the branching problem; fortunately, we do not have that problem in this study. 

Figure 12 shows the GUI (graphic user interface) of this system. The semi-automated and fully 
automated algorithms and the manual tracing function are embedded together. The GUI has zoom 
functions; it thus allows the possibility of manual tracing the sub-pixel resolution to trace the 
boundary. The semi-automated algorithm lets the user select an arbitrary vessel (using ROI selection: 
region of interest) for boundary detection. The software is designed for detecting objects with a 
brighter area than its surroundings. Detecting the dark area is possible by inverting the gradient 
images. 

Figure 12. The software system GUI. The system can read DICOM and other image formats 
supported by MatLab. However, only the DICOM format offers pixel size information. 
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The applications of the proposed algorithm are not restricted to vessel segmentation. It is also 
practical for bone segmentation in 3D CT slices [26] and detecting the boundaries of intima and 
adventitia of common carotid artery inner walls in dynamic sonographic image sequences [6], which 
requires high segmentation accuracy. The advantage is that this method provides the connectivity 
information in three dimensions. Many segmentation methods do not consider the third dimension. The 
computation is also fast in this method, so real-time implementation on the object tracking problem 
might be possible. 

Our future aim is to test this system on sonographic image sequences to detect the intima on both 
sides of the CCA inner wall to extract the inner wall lumen diameter changing with time. Combined 
with the blood pressure information, the lumen diameter change can compute the arterial elasticity, 
which is an important marker in clinics. As a final remark, we must note that this method might be 
extended to the dual surface detection in a 3D volume, such as the dual dynamic programming in a 2D 
case [8]. 
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6. Conclusions 

We have developed a novel algorithm that can detect the optimal surface in a 3D volume matrix. 
This algorithm is applied to detect SFA boundaries in MRA image sequences. According to our 
phantom study, we conclude that the optimal parameter values for using the 3D DP are thus: image 
resize factor = 1.6, and weighting factor s = 0.02. In the MRA real image study, the average relative 
unsigned error is 2.4% ± 2.0% in 16 MRI sequences (800 images) compared to the manual tracings. 
The algorithm proposed in this study is reliable with repeatable results that can replace the experts' 
manual work. 
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